マウス15.5日胚の眼球から取り出した網膜の細胞。10X Chromium 3′ v2で解析。2つのReplicatesでそれぞれ約2,600細胞。
Mouse developing retina from Lo Giudice et al. https://doi.org/10.1242/dev.178103
Lo Giudice, Q., Leleu, M., La Manno, G. & Fabre, P. J. Single-cell transcriptional logic of cell-fate specification and axon guidance in early-born retinal neurons. Development 146, dev178103 (2019).
好きな方法で公共データベースからダウンロードする。レプリケイトがふたつあって、それぞれ以下のSRR IDに対応。
fastq-dump SRR8181428 --split-files --gzip
fastq-dump SRR8181429 --split-files --gzip
Cell Rangerでマウスゲノム(mm10)へのマッピング、遺伝子ごとのカウントを実行。
レプリケイトは(論文でそう名付けられてたので)E2, F2という名前のバッチに設定。
cellranger count --id=RetinalBatchE2 \
--fastqs=/path/to/SRR8181428\
--sample=SRR8181428 \
--transcriptome=/path/to/Cell_Ranger_References/refdata-gex-mm10-2020-A
cellranger count --id=RetinalBatchF2 \
--fastqs=/path/to/SRR8181429\
--sample=SRR8181429 \
--transcriptome=/path/to/Cell_Ranger_References/refdata-gex-mm10-2020-A
講習では、このCell Ranger解析の結果ディレクトリ(の一部、講習で使う部分だけを抜き出したもの)を配布している。
こんな感じの構成になっていて、Cell Rangerのアウトプットが見えるはず。velocytoは後半で解説。
!tree ./data
./data ├── RetinalBatchE2 │ ├── outs │ │ ├── filtered_feature_bc_matrix │ │ │ ├── barcodes.tsv.gz │ │ │ ├── features.tsv.gz │ │ │ └── matrix.mtx.gz │ │ └── web_summary.html │ └── velocyto │ └── RetinalBatchE2.loom ├── RetinalBatchF2 │ ├── outs │ │ ├── filtered_feature_bc_matrix │ │ │ ├── barcodes.tsv.gz │ │ │ ├── features.tsv.gz │ │ │ └── matrix.mtx.gz │ │ └── web_summary.html │ └── velocyto │ └── RetinalBatchF2.loom ├── retinal.h5ad └── retinal_velo.loom 8 directories, 12 files
定番のライブラリのimportに加えて、scanpyをscの別名でimportする。
Figのちょっとした設定。あと一応、筆者の環境における各ライブラリのバージョンは以下の通り。
import warnings
warnings.filterwarnings('ignore')
import os
import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt
import seaborn as sns
sc.logging.print_header()
sc.settings.set_figure_params(dpi=100, facecolor='white')
scanpy==1.9.1 anndata==0.8.0 umap==0.5.3 numpy==1.22.3 scipy==1.9.3 pandas==1.5.2 scikit-learn==1.1.3 statsmodels==0.13.5 python-igraph==0.10.2 pynndescent==0.5.8
Scanpyの特徴は、データフレームを拡張した *anndata* と呼ばれる オブジェクトを使う点にある。anndataを使うと、ひとつのオブジェクトに遺伝子発現量のデータ、サンプルや細胞のアノテーション、遺伝子の情報などをまとめて格納できる。 anndataを使うことで、実験の情報が詰まったひとつのオブジェクトに対して処理を次々に実行し、さらに処理結果をそこに蓄積していくことができる。
Scanpyでは、10x Genomicsのデータは、結果のディレクトリを指定することでそのままロードが可能な関数が用意されている。
ディレクトリには、3つのファイルが書き出されている。
ひとつは、個別の細胞を識別するバーコードが記述されたbarcodes.tsvという名前のテキストファイル。各行ごとにひとつのバーコードが記述されている。
ふたつめは、計測された遺伝子が記述されたgenes.tsvという名前のタブ区切りテキストファイル。このファイルは左側の列に遺伝子のENSEMBL Gene ID、右側の列の遺伝子のシンボルが記述されている。
最後にmatrix.mtxというテキストファイルが、各細胞(バーコード)について各遺伝子の発現がいくつ観測されたのかカウント情報をまとめたファイル。このファイルはMatrix Market Exchange Formatsという形式で記述されており、疎行列(含まれる値の多くがゼロであるような行列)を比較的コンパクトに記述するための形式となっている。
Scanpyの10xデータ用読み込み関数を使うと、これらを同時に読み込んで、適切に構造化されたオブジェクトができあがる。
adata_E2 = sc.read_10x_mtx(path='./data/RetinalBatchE2/outs/filtered_feature_bc_matrix/', cache=True)
adata_F2 = sc.read_10x_mtx(path='./data/RetinalBatchF2/outs/filtered_feature_bc_matrix/', cache=True)
adata_E2
AnnData object with n_obs × n_vars = 3166 × 32285
var: 'gene_ids', 'feature_types'
adata_F2
AnnData object with n_obs × n_vars = 3299 × 32285
var: 'gene_ids', 'feature_types'
2つのバッチをひとつのanndataオブジェクトに統合する。バッチのラベルもここで設定。
adata = adata_E2.concatenate(adata_F2, batch_categories=['E2', 'F2'])
観測値(細胞)に関するデータは、obs でアクセスする。 observationsの略。このデータフレームに細胞に関するデータを追加していくことができる。
adata.obs
| batch | |
|---|---|
| AAACCTGAGATGTCGG-1-E2 | E2 |
| AAACCTGCAATCCAAC-1-E2 | E2 |
| AAACCTGCACATGGGA-1-E2 | E2 |
| AAACCTGCAGCAGTTT-1-E2 | E2 |
| AAACCTGGTTCCTCCA-1-E2 | E2 |
| ... | ... |
| TTTGTCACATCGATGT-1-F2 | F2 |
| TTTGTCAGTAGAGGAA-1-F2 | F2 |
| TTTGTCAGTCATCCCT-1-F2 | F2 |
| TTTGTCATCAGTTCGA-1-F2 | F2 |
| TTTGTCATCCGAATGT-1-F2 | F2 |
6465 rows × 1 columns
変数(遺伝子)に関するデータは、var でアクセスする。variable の略。現在は、遺伝子のシンボル(名前)と、Ensembl Gene IDなどが登録されている。
adata.var
| gene_ids | feature_types | |
|---|---|---|
| Xkr4 | ENSMUSG00000051951 | Gene Expression |
| Gm1992 | ENSMUSG00000089699 | Gene Expression |
| Gm19938 | ENSMUSG00000102331 | Gene Expression |
| Gm37381 | ENSMUSG00000102343 | Gene Expression |
| Rp1 | ENSMUSG00000025900 | Gene Expression |
| ... | ... | ... |
| AC124606.1 | ENSMUSG00000095523 | Gene Expression |
| AC133095.2 | ENSMUSG00000095475 | Gene Expression |
| AC133095.1 | ENSMUSG00000094855 | Gene Expression |
| AC234645.1 | ENSMUSG00000095019 | Gene Expression |
| AC149090.1 | ENSMUSG00000095041 | Gene Expression |
32285 rows × 2 columns
カウントテーブルには、以下のように *X* でアクセスできる。疎行列として格納されている。
adata.X
<6465x32285 sparse matrix of type '<class 'numpy.float32'>' with 10273670 stored elements in Compressed Sparse Row format>
adata
AnnData object with n_obs × n_vars = 6465 × 32285
obs: 'batch'
var: 'gene_ids', 'feature_types'
観測された遺伝子が極端に少ない細胞、割り当てられた細胞が極端に少ない遺伝子をフィルタリング。
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=20)
adata
AnnData object with n_obs × n_vars = 6462 × 13332
obs: 'batch', 'n_genes'
var: 'gene_ids', 'feature_types', 'n_cells'
ミトコンドリア遺伝子の発現割合でフィルタリングする予定なので、遺伝子側のメタデータにミトコンドリア遺伝子か否かの情報を付与。
adata.var['mt'] = adata.var_names.str.startswith('mt-')
以下は、クオリティコントロールに関連したいくつかの指標を一挙に計算してくれる便利な関数。
qc_varsに遺伝子側のメタデータラベルを設定すると、それについてもカウントとパーセンテージを勝手に計算してくれる。
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
adata.obs
| batch | n_genes | n_genes_by_counts | total_counts | total_counts_mt | pct_counts_mt | |
|---|---|---|---|---|---|---|
| AAACCTGAGATGTCGG-1-E2 | E2 | 1432 | 1427 | 3263.0 | 111.0 | 3.401778 |
| AAACCTGCAATCCAAC-1-E2 | E2 | 1297 | 1292 | 2568.0 | 97.0 | 3.777258 |
| AAACCTGCACATGGGA-1-E2 | E2 | 261 | 261 | 572.0 | 47.0 | 8.216784 |
| AAACCTGCAGCAGTTT-1-E2 | E2 | 360 | 359 | 624.0 | 21.0 | 3.365385 |
| AAACCTGGTTCCTCCA-1-E2 | E2 | 3072 | 3057 | 10505.0 | 283.0 | 2.693955 |
| ... | ... | ... | ... | ... | ... | ... |
| TTTGTCACATCGATGT-1-F2 | F2 | 1916 | 1908 | 3714.0 | 136.0 | 3.661820 |
| TTTGTCAGTAGAGGAA-1-F2 | F2 | 2010 | 2007 | 5526.0 | 176.0 | 3.184944 |
| TTTGTCAGTCATCCCT-1-F2 | F2 | 1760 | 1757 | 3718.0 | 121.0 | 3.254438 |
| TTTGTCATCAGTTCGA-1-F2 | F2 | 1175 | 1170 | 2064.0 | 54.0 | 2.616279 |
| TTTGTCATCCGAATGT-1-F2 | F2 | 443 | 441 | 865.0 | 305.0 | 35.260117 |
6462 rows × 6 columns
それぞれの分布をプロット。
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
jitter=0.4, multi_panel=True)
sc.pl.scatter(adata, 'total_counts', 'n_genes_by_counts', color='pct_counts_mt', size=40)
fig = plt.figure()
sns.displot(adata.obs['pct_counts_mt'][adata.obs['pct_counts_mt'] < 10], kde=False)
plt.show()
<Figure size 400x400 with 0 Axes>
分布を見て方針を決めて、フィルタリングを実行。
print('Total number of cells: {:d}'.format(adata.n_obs))
sc.pp.filter_cells(adata, min_counts = 2000)
print('Number of cells after min count filter: {:d}'.format(adata.n_obs))
sc.pp.filter_cells(adata, max_counts = 13000)
print('Number of cells after max count filter: {:d}'.format(adata.n_obs))
sc.pp.filter_cells(adata, min_genes = 1000)
print('Number of cells after gene filter: {:d}'.format(adata.n_obs))
adata = adata[adata.obs['pct_counts_mt'] < 6]
print('Number of cells after MT filter: {:d}'.format(adata.n_obs))
Total number of cells: 6462 Number of cells after min count filter: 4728 Number of cells after max count filter: 4706 Number of cells after gene filter: 4686 Number of cells after MT filter: 4563
adata
View of AnnData object with n_obs × n_vars = 4563 × 13332
obs: 'batch', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'n_counts'
var: 'gene_ids', 'feature_types', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
正規化や対数変換などの前に、この段階のカウントマトリックスを別のレイヤーに退避させておく。あとで確率モデル推定のときに必要になる。レイヤーに保管したデータは基本的に以降の操作の影響を受けないが、細胞や遺伝子のslicingは適用されるので注意。(常にマトリックスのshapeは一致する)
adata.layers['counts'] = adata.X.copy()
ライブラリサイズによる正規化、対数変換などの前処理は、scanpy.pp以下にいくつか便利な関数がある。
ここでは、細胞ごとのカウントの和が10,000になるように正規化してから、対数変換。
sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)
sc.pp.log1p(adata)
この段階のデータ(いろいろとややこしい変換をしていない「ピュア」なデータ)も、あとでプロットや結果解釈のときに使いたいので退避させておく。rawは特別なレイヤーで、数値だけでなくobsやvarのメタデータも含めてフリーズしてくれる。また遺伝子側のsliceの影響を受けない(細胞側のsliceは適用される)
adata.raw = adata
発現量変動の大きい遺伝子のみを抽出して、データのサイズを小さくする。
sc.pp.highly_variable_genes(adata, flavor='seurat_v3', n_top_genes=2000)
print('\n','Number of highly variable genes: {:d}'.format(np.sum(adata.var['highly_variable'])))
Number of highly variable genes: 2000
sc.pl.highly_variable_genes(adata)
計算結果は自動的に adata.var 、つまり、遺伝子に関するメタデータを格納したオブジェクトに追加される。 highly_variable の項目が True の遺伝子が高発現変動遺伝子。
adata.var
| gene_ids | feature_types | n_cells | mt | n_cells_by_counts | mean_counts | pct_dropout_by_counts | total_counts | highly_variable | highly_variable_rank | means | variances | variances_norm | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Xkr4 | ENSMUSG00000051951 | Gene Expression | 89 | False | 89 | 0.013773 | 98.622717 | 89.0 | False | NaN | 0.013526 | 0.013121 | 0.834873 |
| Gm19938 | ENSMUSG00000102331 | Gene Expression | 578 | False | 578 | 0.103528 | 91.055401 | 669.0 | False | NaN | 0.118130 | 0.129846 | 0.984716 |
| Mrpl15 | ENSMUSG00000033845 | Gene Expression | 1899 | False | 1899 | 0.375271 | 70.612813 | 2425.0 | False | NaN | 0.471746 | 0.418216 | 0.943994 |
| Lypla1 | ENSMUSG00000025903 | Gene Expression | 1284 | False | 1284 | 0.241566 | 80.129991 | 1561.0 | False | NaN | 0.302353 | 0.307116 | 0.989431 |
| Tcea1 | ENSMUSG00000033813 | Gene Expression | 2614 | False | 2614 | 0.631848 | 59.548128 | 4083.0 | False | NaN | 0.704444 | 0.549619 | 0.948720 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| CAAA01118383.1 | ENSMUSG00000063897 | Gene Expression | 96 | False | 96 | 0.015166 | 98.514392 | 98.0 | False | NaN | 0.021325 | 0.025196 | 1.020294 |
| Vamp7 | ENSMUSG00000051412 | Gene Expression | 1606 | False | 1606 | 0.314144 | 75.147013 | 2030.0 | False | NaN | 0.395332 | 0.365986 | 0.947294 |
| Tmlhe | ENSMUSG00000079834 | Gene Expression | 46 | False | 46 | 0.007119 | 99.288146 | 46.0 | False | NaN | 0.007515 | 0.008047 | 0.921716 |
| 4933409K07Rik | ENSMUSG00000095552 | Gene Expression | 20 | False | 20 | 0.003095 | 99.690498 | 20.0 | False | NaN | 0.003088 | 0.002459 | 0.797902 |
| AC149090.1 | ENSMUSG00000095041 | Gene Expression | 1960 | False | 1960 | 0.475085 | 69.668833 | 3070.0 | False | NaN | 0.462647 | 0.455644 | 1.044050 |
13332 rows × 13 columns
PCAはscanpyの前処理関数で簡単に実行できる。とりあえず、50次元まで落としてみる。 use_highly_variableのフラグをオンにすると、遺伝子全体ではなく前項で決定した高発現変動遺伝子のみを使って次元削減をする。
sc.pp.pca(adata, n_comps=50, use_highly_variable=True, svd_solver='arpack')
PCAで落とした座標は、観察値のメタデータを格納する *obsm* に自動的に入る。obsm は Multi-dimensional annotation of observations の略。複数の次元でひとつの何かを表現するような観測値のメタデータがここに格納される。
print(adata.obsm['X_pca'].shape)
(4563, 50)
プロットしてみる。Scanpyでは基本的に「前処理」(Preprocessing)に関わる関数がscanpy.ppに、「プロット」(Plot)に関わる関数がscanpy.plに入っている。
sc.pl.pca(adata, color='total_counts')
sc.pl.pca(adata, color='total_counts',
components=['1,2', '2,3', '1,3'])
まず、t-SNE, UMAP共通のステップとして、データから「近傍グラフ」(neighborhood graph)の構築が必要。
sc.pp.neighbors(adata)
データ点間の接続関係(細胞の近傍関係)は、全細胞vs.全細胞のペアの情報を記録する *obsp* (Pairwise annotation of observations の略)に格納される。
adata.obsp['connectivities']
<4563x4563 sparse matrix of type '<class 'numpy.float32'>' with 94424 stored elements in Compressed Sparse Row format>
接続関係を元にしてt-SNEを実行。
sc.tl.tsne(adata)
t-SNE結果のプロット。いまのところ特に色つけるようなデータも無いので、適当に細胞ごとのカウントで色分けをした。遺伝子名をcolorに指定するとその遺伝子の発現量によるプロットも可能。
sc.pl.tsne(adata, color='total_counts')
近傍グラフはすでに計算しているので、scanpy.tlのumap関数を使えばオーケー。
sc.tl.umap(adata)
プロットはt-SNEとほとんど同じ。scanpy.pl以下にUMAP用の描画関数がある。
sc.pl.umap(adata, color='total_counts')
あきらかに似た形状の、ふたつの構造がある。。。
各細胞が由来するバッチで色分けしてみると・・・
sc.pl.umap(adata, color='batch')
ということで、このふたつの構造はもろにバッチ効果の結果だった。なんらかの方法で補正が必要だが、それは後述。
Leidenクラスタリングを実行してみる。モジュラリティの計算に影響を与える "resolution" パラメータが存在する。
sc.tl.leiden(adata, resolution=0.5, key_added='leiden_r0.5')
クラスタリングの結果は観測値のメタデータ(obs)の中に格納される。
adata.obs
| batch | n_genes | n_genes_by_counts | total_counts | total_counts_mt | pct_counts_mt | n_counts | leiden_r0.5 | |
|---|---|---|---|---|---|---|---|---|
| AAACCTGAGATGTCGG-1-E2 | E2 | 1427 | 1427 | 3263.0 | 111.0 | 3.401778 | 3263.0 | 1 |
| AAACCTGCAATCCAAC-1-E2 | E2 | 1292 | 1292 | 2568.0 | 97.0 | 3.777258 | 2568.0 | 1 |
| AAACCTGGTTCCTCCA-1-E2 | E2 | 3057 | 3057 | 10505.0 | 283.0 | 2.693955 | 10505.0 | 2 |
| AAACCTGTCCAATGGT-1-E2 | E2 | 1449 | 1449 | 2781.0 | 91.0 | 3.272204 | 2781.0 | 7 |
| AAACGGGAGGCAATTA-1-E2 | E2 | 2772 | 2772 | 9240.0 | 224.0 | 2.424242 | 9240.0 | 1 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| TTTGTCAAGGGCTCTC-1-F2 | F2 | 2393 | 2393 | 6209.0 | 143.0 | 2.303108 | 6209.0 | 5 |
| TTTGTCACATCGATGT-1-F2 | F2 | 1908 | 1908 | 3714.0 | 136.0 | 3.661820 | 3714.0 | 6 |
| TTTGTCAGTAGAGGAA-1-F2 | F2 | 2007 | 2007 | 5526.0 | 176.0 | 3.184944 | 5526.0 | 0 |
| TTTGTCAGTCATCCCT-1-F2 | F2 | 1757 | 1757 | 3718.0 | 121.0 | 3.254438 | 3718.0 | 5 |
| TTTGTCATCAGTTCGA-1-F2 | F2 | 1170 | 1170 | 2064.0 | 54.0 | 2.616279 | 2064.0 | 6 |
4563 rows × 8 columns
プロットする点の座標じたいはUMAPの結果なので、プロットはUMAP版の関数を使い、色分けだけをクラスタリング結果で指定すればいい。
sc.pl.umap(adata,
color=['batch', 'leiden_r0.5'],
ncols=2,
frameon=False)
scVIが仮定しているデータの生成プロセスは以下の確率モデルに基づく。
全細胞数がN, 細胞のインデックスが $n \in \{1..N\}$、全遺伝子数がG、遺伝子のインデックスが $g \in \{1..G\}$、細胞nのバッチIDが$s_n$であるとしたとき、細胞nの遺伝子gのカウントは、
$$z_n \sim Normal(0, I)$$$$l_n \sim logNormal(l_\mu, l_\sigma^2)$$$$\rho_n = f_w (z_n, s_n)$$$$w_{ng} = Gamma(\rho_n^g, \theta)$$$$y_{ng} = Poisson(l_n w_{ng})$$$$h_{ng} = Bernoulli( f_h^g(z_n, s_n) )$$$$ \begin{equation} x_{ng} = \begin{cases} y_{ng} & \text{if}\ h_{ng}=0 \\ 0 & \text{otherwise} \end{cases} \end{equation} $$$l_\mu, l_\sigma$は、スケーリングファクターでバッチの数をBとしたとき$l_\mu,l_\sigma \in \mathbb{R}^B$で、実際のバッチごとの平均と分散に設定しておく。
ガンマ分布のパラメータ、ベルヌーイ分布のパラメータは、正規分布からサンプリングした$z_n$を使ってニューラルネットワークで表現する(VAEにおけるreparametrization trick)。
最終的なカウントはガンマとポアソンの複合分布なので、負の二項分布でモデル化してることになる。さらにベルヌーイ分布で、scRNA-seqにありがちなdrop-outイベントを表現。これらの組み合わせで、ゼロ過剰負の二項分布(Zero-inflated negative binomial distribution; ZINB)を表現している。
事後分布はVAEをSGDで最適化して変分推論。
学習後、$z_n$にアクセスすれば、バッチの効果が除去された細胞の潜在ベクトルが得られるし、$\rho_{ng}$にアクセスすればdrop-outなどのイベントが生じなかった場合の遺伝子発現の期待値が得られる。
import scvi
Global seed set to 0 /root/anaconda3/envs/pags/lib/python3.9/site-packages/pytorch_lightning/utilities/warnings.py:53: LightningDeprecationWarning: pytorch_lightning.utilities.warnings.rank_zero_deprecation has been deprecated in v1.6 and will be removed in v1.8. Use the equivalent function from the pytorch_lightning.utilities.rank_zero module instead. new_rank_zero_deprecation( /root/anaconda3/envs/pags/lib/python3.9/site-packages/pytorch_lightning/utilities/warnings.py:58: LightningDeprecationWarning: The `pytorch_lightning.loggers.base.rank_zero_experiment` is deprecated in v1.7 and will be removed in v1.9. Please use `pytorch_lightning.loggers.logger.rank_zero_experiment` instead. return new_rank_zero_deprecation(*args, **kwargs)
データを格納したanndataオブジェクトと、バッチのIDを指定したobsのカラムのラベルを指定して、モデルをセットアップする。カウントデータの確率モデルなので、(正規化や対数変換したマトリックスではなく)カウントデータを格納したレイヤーを指定する必要があることに注意。
scvi.model.SCVI.setup_anndata(
adata,
layer='counts',
batch_key='batch',
)
WARNING:jax._src.lib.xla_bridge:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
モデルを作る。ここで中間層のノード数とか潜在ベクトルの次元サイズなど設定できるが、とくに変更しなくていいと思う。
model = scvi.model.SCVI(adata)
学習を実行する。GPU使わない場合はめっちゃ時間かかるので今回の講習では割愛。学習済みのモデルパラメータを用意したのでそれをロードして学習できたことにする。
# model.train()
# model.save('./models/scVI_model', overwrite=True)
model = scvi.model.SCVI.load('./models/scVI_model', adata=adata)
INFO File ./models/scVI_model/model.pt already downloaded
学習された細胞ごとの潜在表現$z_n$は以下の関数で取得できる。PCAやUMAPの座標と同じように、obsmに格納しておく。
adata.obsm['X_scVI'] = model.get_latent_representation()
発現量期待値 $\rho_{ng}$ 、つまりdenoiseされた発現量テーブルは以下の関数で取得できる。適当なライブラリサイズで規格化してくれる。この数値はあとで使いたいので別のレイヤーに格納しておく。
adata.layers['scvi_normalized'] = model.get_normalized_expression(library_size=1e4)
adata
AnnData object with n_obs × n_vars = 4563 × 13332
obs: 'batch', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'n_counts', 'leiden_r0.5', '_scvi_batch', '_scvi_labels'
var: 'gene_ids', 'feature_types', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm'
uns: 'log1p', 'hvg', 'pca', 'neighbors', 'tsne', 'umap', 'batch_colors', 'leiden', 'leiden_r0.5_colors', '_scvi_uuid', '_scvi_manager_uuid'
obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_scVI'
varm: 'PCs'
layers: 'counts', 'scvi_normalized'
obsp: 'distances', 'connectivities'
潜在表現でちゃんとバッチ効果が補正されたかどうか確認してみる。
近傍グラフ計算の scanpy.pp.neighbors は、デフォルトではPCAで計算した主成分(adata.obsm['X_pca'])を元にグラフを構築する。ここではPCA結果の座標ではなく、scVIで推定された細胞の潜在ベクトルを指定してグラフ構築してみる。バッチ効果が補正されたグラフが構築されるはず。
sc.pp.neighbors(adata,
n_neighbors=30,
use_rep="X_scVI")
sc.tl.umap(adata, min_dist=0.5)
sc.pl.umap(adata, color='batch')
同じように、scVI潜在表現で構成された近傍グラフを元にLeidenクラスタリングを実行してみる。上のセルで scanpy.pp.neighbors を実行したことにより、adata.obspにはscVI潜在表現で構築された近傍グラフが入っているので、ここでは普通に scanpy.tl.leiden を実行すればいい。
sc.tl.leiden(adata, key_added="leiden_scVI", resolution=1.0)
sc.pl.umap(adata,
color=['batch', 'leiden_scVI'],
ncols=2,
frameon=False)
SOLOによるDoublet検出。
Bernstein, Nicholas J., et al. "Solo: doublet identification in single-cell RNA-Seq via semi-supervised deep learning." Cell systems 11.1 (2020): 95-101.
時間に余裕があればやる。
results = []
for batch in ['F2', 'E2']:
tmp_solo_model = scvi.external.SOLO.from_scvi_model(
model, restrict_to_batch=batch)
tmp_solo_model.train()
result = tmp_solo_model.predict(soft=False)
# 意図は不明だがSOLO予測のindexになぜか "-0" が付加されるので消しておく
result.index = result.index.str.replace("-0$", "", regex=True)
results.append(result)
results = pd.concat(results)
INFO Creating doublets, preparing SOLO model.
GPU available: False, used: False TPU available: False, using: 0 TPU cores IPU available: False, using: 0 IPUs HPU available: False, using: 0 HPUs
Epoch 79/400: 20%|█▉ | 79/400 [00:13<00:54, 5.90it/s, loss=0.144, v_num=1]
Monitored metric validation_loss did not improve in the last 30 records. Best score: 0.176. Signaling Trainer to stop.
INFO Creating doublets, preparing SOLO model.
GPU available: False, used: False TPU available: False, using: 0 TPU cores IPU available: False, using: 0 IPUs HPU available: False, using: 0 HPUs
Epoch 104/400: 26%|██▌ | 104/400 [00:15<00:45, 6.56it/s, loss=0.131, v_num=1] Monitored metric validation_loss did not improve in the last 30 records. Best score: 0.158. Signaling Trainer to stop.
adata.obs['SOLO_prediction'] = results[adata.obs.index]
adata.obs['SOLO_prediction'].value_counts()
singlet 4109 doublet 454 Name: SOLO_prediction, dtype: int64
sc.pl.umap(adata, color='SOLO_prediction', frameon=False)
仮説検定ではなく、ベイズ統計の枠組みでモデル比較を行う。基本的にはある集団と別の集団それぞれの遺伝子発現量期待値(VAEでdenoiseした発現量$\rho_{ng}$)のlog fold changeの値が、ある閾値以下となるか、それ以上の変化があるか、のふたつのモデルで比較する。
groupbyにラベルを指定すると、そのクラス(この場合はleidenクラスタリングのひとつのクラスタ)とその他全部で自動的に比較してくれる。比較のグループを一部に制限したり、ペアを指定したりいろいろと複雑に組み合わせは指定できる。詳細はAPIのドキュメントに。
DEs = model.differential_expression(groupby='leiden_scVI')
DEs.head()
DE...: 100%|██████████| 11/11 [00:54<00:00, 4.95s/it]
| proba_de | proba_not_de | bayes_factor | scale1 | scale2 | pseudocounts | delta | lfc_mean | lfc_median | lfc_std | ... | raw_mean1 | raw_mean2 | non_zeros_proportion1 | non_zeros_proportion2 | raw_normalized_mean1 | raw_normalized_mean2 | is_de_fdr_0.05 | comparison | group1 | group2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Carmil2 | 0.9836 | 0.0164 | 4.093937 | 2.658159e-07 | 0.000005 | 0.0 | 0.25 | -5.505758 | -5.591664 | 5.274880 | ... | 0.001264 | 0.015907 | 0.001264 | 0.015642 | 0.001908 | 0.021632 | True | 0 vs Rest | 0 | Rest |
| Pcdhb18 | 0.9822 | 0.0178 | 4.010596 | 8.968271e-08 | 0.000006 | 0.0 | 0.25 | -6.042808 | -6.083243 | 6.751517 | ... | 0.000000 | 0.006628 | 0.000000 | 0.006098 | 0.000000 | 0.012190 | True | 0 vs Rest | 0 | Rest |
| Cntn5 | 0.9812 | 0.0188 | 3.954919 | 2.587765e-07 | 0.000017 | 0.0 | 0.25 | -5.148154 | -5.396956 | 4.268237 | ... | 0.001264 | 0.074231 | 0.001264 | 0.056469 | 0.002287 | 0.091536 | True | 0 vs Rest | 0 | Rest |
| Gfra3 | 0.9808 | 0.0192 | 3.933458 | 2.856224e-08 | 0.000010 | 0.0 | 0.25 | -6.871420 | -7.264495 | 5.107664 | ... | 0.000000 | 0.022269 | 0.000000 | 0.016702 | 0.000000 | 0.037028 | True | 0 vs Rest | 0 | Rest |
| Psd2 | 0.9808 | 0.0192 | 3.933458 | 9.740267e-08 | 0.000005 | 0.0 | 0.25 | -6.067220 | -6.589231 | 5.682239 | ... | 0.000000 | 0.016437 | 0.000000 | 0.015642 | 0.000000 | 0.023997 | True | 0 vs Rest | 0 | Rest |
5 rows × 22 columns
それぞれのクラスタ「特異的」遺伝子、DEG確率の高い順に上から3つくらい見てみる。
markers = {}
for i, c in enumerate(adata.obs['leiden_scVI'].unique()):
comparison_label = "{} vs Rest".format(c)
cluster_df = DEs.loc[DEs['comparison'] == comparison_label]
cluster_df = cluster_df[cluster_df['lfc_mean'] > 0]
cluster_df = cluster_df[cluster_df['bayes_factor'] > 3]
cluster_df = cluster_df[cluster_df['non_zeros_proportion1'] > 0.1]
markers[c] = cluster_df.index.tolist()[:3]
図でLeidenクラスタのデンドログラムを表示したいので事前に計算しておく。
sc.tl.dendrogram(adata, groupby='leiden_scVI', use_rep='X_scVI')
クラスタ間での遺伝子発現量の比較は、scanpy.pl.dotplotが見やすくて便利。上で抽出したマーカー遺伝子を指定して、表示する発現量の数値としてはadata.rawに格納したデータで描画する。
sc.pl.dotplot(
adata,
markers,
groupby='leiden_scVI',
dendrogram=True,
color_map="Blues",
swap_axes=True,
use_raw=True,
standard_scale='var')
クラスタ平均値ではなくクラスタごと細胞ごとの全部の(scVIノーマライズされた)発現量でヒートマップを描くこともできる。
sc.pl.heatmap(
adata,
markers,
groupby='leiden_scVI',
layer='scvi_normalized',
standard_scale="var",
dendrogram=True,
figsize=(8, 12)
)
adata.write(filename='./data/retinal.h5ad')